; docformat = 'rst'
;
;+
;
; :Purpose:
;   Calculates emissions within each LPDM region, based on global emissions file
;   {input_directory}/emissions/emissions.inp.sav
;   
; :Inputs:
;   Sim_name: (string) Simulation name
;
; :Requires:
;   Martin Shultz's REGRID routines. Download at http://www.acd.ucar.edu/Applications/regrid_v2.tar.gz
;   Matt Rigby's IDL library: http://code.google.com/p/mr-idl-library/
;
; :Outputs:
;   Generates an emissions file for each region:
;   
;   E.g. {input_directory}/{sim_name}/emissions/lr_00_08.sav contains emissions for the first (index 00) local region in 2008
;   
; :History:
; 	Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
pro cels_emissions_lpdm, sim_name

  compile_opt idl2
  on_error, 2
  
  print, "CELS_EMISSIONS_LPDM: LPDM emissions"
  
  stations=cels_get_parameter(sim_name, 'STATIONS')
  StartY=cels_get_parameter(sim_name, 'STARTY')
  EndY=cels_get_parameter(sim_name, 'ENDY')

  nStations=n_elements(stations)

  restore, cels_filestr(/Input, sim_name + '/emissions/emissions.inp.sav')
  restore, cels_filestr(/Input, sim_name + '/regions.sav')

  ;Make sure emissions grid is set up correctly
  mr_grid_orient, q_lon, q_lat, q

  ;This needs changing if timesteps other than years are used
  caldat, q_time, dummy, dummy1, year
  
  q_input=temporary(q)
  
  for yi=StartY, EndY-1 do begin
  
    wh=where(year eq yi, count)
    if count eq 0 then begin
      message, strcompress('No emissions for year ' + string(yi))
    endif else begin
      q_annual=reform(q_input[*, *, wh])
    endelse
  
    for ri=0, N_LR-1 do begin
      
      lon=reform(lon_local[ri])
      lat=reform(lat_local[ri])

      ;Martin Shultz's regridder
      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
      
      ;Define old grid
      nlon = n_elements(q_lon)
      incr = (max(q_lon) - min(q_lon))/Double(nlon-1)
      ms_lon = Obj_New('MGS_LonVector', $
                    start=min(q_lon),incr=incr,nval=nlon, $
                    type='regular')
      nlat = n_elements(q_lat)
      incr = (max(q_lat) - min(q_lat))/Double(nlat-1)
      ms_lat = Obj_New('MGS_LatVector', $
                    start=min(q_lat),incr=incr,nval=nlat, $
                    type='regular')
      o = Obj_New('MGS_RGrid', name='Source Grid', $
                  Longitude=ms_lon, Latitude=ms_lat)
      if total(o->getlon() - q_lon) gt 0.1 then stop
      if total(o->getlat() - q_lat) gt 0.1 then stop

      ;Define new grid
      nlon = n_elements(lon)
      lon_incr = (max(lon) - min(lon))/Double(nlon-1)
      ms_lon = Obj_New('MGS_LonVector', $
                    start=min(lon),incr=lon_incr,nval=nlon, $
                    type='regular')
      nlat = n_elements(lat)
      lat_incr = (max(lat) - min(lat))/Double(nlat-1)
      ms_lat = Obj_New('MGS_LatVector', $
                    start=min(lat),incr=lat_incr,nval=nlat, $
                    type='regular')
      n = Obj_New('MGS_RGrid', name='Target Grid', $
                  Longitude=ms_lon, Latitude=ms_lat)

      if total(n->getlon() - lon) gt 0.1 then stop
      if total(n->getlat() - lat) gt 0.1 then stop

      ;Regrid
      q=regrid(q_annual, oldgrid=o, newgrid=n, /use_grid, /area_weighted) ;kg/m2/s

      ;Clean
      Obj_Destroy, o
      Obj_Destroy, n

      ;Save output
      lonq=lon
      latq=lat
      save, filename=cels_filestr(/Input, sim_name + '/emissions/lr_' + string(ri, format='(I02)') + '_' + string(yi) + '.sav'), $
        lonq, latq, q

    endfor

  endfor

end